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1.0 INTRODUCTION 


The objective of the work described in this report was the development of a 
finite element scheme applicable to the acoustics of aero-engine ducts. The 
finite element method, which is widely used with considerable success in 
structural mechanics, is currently finding applications in all areas of 
continuum mechanics. 


The application of finite element methods to acoustics is not new. A signifi- 
cant amount of work has already been performed in determining the acoustic 
modes of closed cavities and in determining the response of coupled acoustic 
structural systems (references 1-9). These analyses assumed simple acous 
tic systems in a uniform medium with no mean flow. Under these circum- 
stances it was possible to formulate variational principles which were then 
used as the basis for a finite element model. None of these assumptions 
are valid within a modern- aircrqfFengine such as shown L in Figure 1. Here, 
the noise generation mechanisms are complex, and the medium which trans- 
mits the sound contains large gradients and moves with a considerable 
velocity. Since no variational principal exists for acoustic fluctuations in 
such a medium, an acoustic analysis must proceed directly from the 
differential equations which describe compressible flow. 

From a cursory examination of the physical size of fan engine ducts on 
current aircraft, and from the knowledge that predominant fan tones exceed 
2000Hz, it is clear that the principal difficulty in development of an acoustic 
model on existing computers is the sheer scale of the required model. In 
structural mechanics after 20 years of experience, 10 000 degrees of 
freedom still represent a large problem. At the outset, it is apparent that 
the acoustic model will require upwards of 100 000 degrees of freedom. 

This is not the only complication, because in structural mechanics as well 




as in previous acoustic studies, global finite element matrices possess 
desirable properties, being generally real, symmetric, and positive 
definite. The symmetric, positive=definite form results from the variational 
formulation on which the finite element models are based. However, the 
Galerkin variation of the method of weighted residuals used for the acoustic 
element in this work produces matrices which are complex, non-Hermitian 
and nonpositive definite. Matrices of this type require special techniques 
to avoid numerical errors in their solution. These techniques are less 
plentiful than for symmetric positive -definite matrices and become difficult 
to implement efficiently for large matrix orders. 

This report describes a solution to the problems outlined above, describes 
the development of a finite element model, and shows typical results 
derived from this model. 



2.0 DESCRIPTION OF MATHEMATICAL MODEL, 


2. 1 Problem Formulation 

The purpose of this analysis is the calculation of sound propagation in 
axisymmetric ducts containing compressible flow. The nondimens ional 
equations governing the motion of the fluid are: 


Conservation of Momentum 

( F.v) ) * -vf +J_v.r 
Mat / R e 

Conservation of Mass 

— * ( p y ) - o 

at 


(i) 


( 2 ) 


Conservation of Energy 


fjir + V.vfj - + V.v)J = J_ |_L v. [kvT) + (r-i) 


(3) 


Equation of State 


= f T 


(4) 


Equations (1) through (4) were made nondimensional for simplicity by using 
an arbitrary characteristic length; the ambient values of temperature, density, 
and speed of sound; and arbitrary characteristic values of viscosity and 
thermal conductivity. These equations govern both the mean flow and imposed 
acoustic motion within the duct. 


4 



While determination of the mean flow is necessary for solution of the acoustic 
problem, the assumptions and techniques used to obtain the mean solution are 
not part of this paper. Instead, it is assumed that the mean flow has already 
been obtained and is in suitable form for inclusion in the acoustic problem. 

At the same time, it is necessary to note that some terms in the equations 
(those involving viscosity and heat conduction) which may be important in 
deriving the mean flow solution are not of primary significance in deter- 
mining the acoustic motion. Thus, while these terms might be retained 
for solution of the mean flow problem, they are eliminated here at the outset 
to simplify derivation of the equations describing the acoustic motion. 

To begin this derivation, we write equation (1) in dimensional component 
form in cylindrical coordinates (ignoring viscosity). 


Axial Direction 
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2U + 
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dr r d 9 d z ) 



( 5 ) 


Radial Direction 
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( 7 ) 
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Equation (2) becomes: 

at + ua_f+ V 3f + wil + p3U + a^ + j.3V) + (7 = 0 

di dz dr de , 2z dr r ddj r ^ 

Assuming that all dependent variables are the sum of two motions; one the 
steady mean flow and the other a fluctuating acoustic motion, we write 

T = f(r,e,2) +^> / (r. e,z,t) 
j> = J (r, e, z) + j'( r < e,z,t) 

M = IA (r, e , z) + u'(r, e,z,t) (9) 

V = V' (r, 9, z) + '''(r , 8 , z,t) 

W = W [r, e , i) + w'(r t e, z,t) 

Substituting equations (9) into equations (5) through (8), we derive two sets of 
equations: the first by assuming that the primed quantities are zero, and the 
second by using the full form of equation (9). Subtracting the first set from 
the second, ignoring second-order terms in the primed quantities and assuming 
mean swirl (W ) to be zero, we get 


f + ($v' + $'V)$_U + f V d \jl 

at dr dr 

t ( Ju't fu) 9U + J u Hi + iif » o 

d 2 dz dz 
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d t 


dr 


dr 


(ll) 


+ ( 5V + s'U) *v + j u **'+ i£ 

d z dz dr 


~ O 


i*' + vi 

/ w'\ 

+ U 2uf' 

st l 

dr r ) 

d Z . 


ss 0 


(12) 


d S + Li d_$ + V d_S + tf/ + V ' B $ 

dt d 2 3r 3z 


+ / 9 n ' + + J_ 9 U2 ' + ✓ x \ 

[ ar r 3e r ) 

+ f'/ sU + 9 V + 7 \ = 0 

l 5 2 5r r / 


The neglect of viscosity and heat conduction in the acoustic motion means 
that it is an adiabatic process. When such is true, the energy equation and 
equation of state may be manipulated to give (Ref. 10, IT) 



(14) 


Here C is the local speed of sound which is determined from the mean flow 


by the relation 



(15) 
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Assuming that C - C (r, z) , we may use equation (14) to eliminate $ -from 
equations (10) through (13). Taking a harmonic solution in t and © for 
the fluctuating quanties 



( u , v, iv, /> ) 


u) i + m 0 

e 


(16) 


where lO and m are pure imaginary, the four equations governing the 
linearized acoustic motion in a nonuniform axisymmetric duct become 


[*.hV ’ 


+ 1/ du + 

I u +_kJL 

| 9 U + U du_ + 

= o 

sW 

dr 

dr 

l U x j 

32 dZ 

$ Bz 


(17) 


Qv +( v 4- 1 Y„ ] U/ + Vdv + [u +1JL) JlY + Udv +j_.Ah - o 

\ S c x ) dr dr \ Sc x J dz fa 3 dr 

(18) 

odw + V I dw + w] -t U <tjv + rn b » o 
l dr r J dz Sr' 


< 2 -b + JLflt ~±th ] +J Lilt - ±t 

c a l dz c dz/ c a ( dr c dr) 

+ v + J / + 2 / + 2*2. I*' + + t I dV + 3 V + / | — q 

dr l dz dr r rj c \ 3z dr r ) 

( 20 ) 




2 .. 2 Solution Strategy 

Solution of equations (17) through (20) is not possible algebraically except 
with severely restrictive assumptions. The question is thus not whether to 
use a numerical technique for their solution, but rather which numerical 
technique to use. 


Within each branch of applied science, partly from its historical develop- 
ment and partly from rational grounds, there appears to exist a favorite 
numerical method. Perhaps it is fortuitous that the rebirth of acoustics is 
too recent to have allowed development of any such prejudice. 


For example, a variety of numerical techniques are used in connection 


with the problem of sound propagation in ducts. Mungur, in his work with 

//•. (• fot* -/•- ('l-t)-^ e 

/Plumblee/\ 1 and .Gladwell » . used a fourth-order Runge-Kutta approach. 

A UHiS- s e%m 

More recently, Baumeister, with Bittner^ s t \ and with Rice^ /^.used an 
integral formulation to generate finite difference approximations to the 
differential equations. A^m o^e indirect technique involving modal expansions 

, . » , ■ ('l 4', 4:&) 

was used by Zorumski 'r A n . . . . , , 

Alsoun this vein is the wave envelope tech- 


( 1 £ ^ 


nique of Nayfeh, Shaker, and Kaiser 
method of multiple scales of the earlier work of Nayfeh and Telionis 


which reduces in the limit to the 

(17) 


Direct numerical solutions are subjept to the problem of size. As pointed 
out by Nayfeh, Shaker, and Kaiser for a negative Mach number 

tending to unity in magnitude, and high frequency, the number of axial steps 
required to resolve the quasi- sinusoidal spatial variation of the fluctuating 
parameters becomes very large. 




No study has yet been conducted, however, to accurately determine the 
limits of direct numerical solutions using currently available computers. 


- 9 /; 



This work began with the objective of determining these limits, using care- 
fully determined storage and matrix manipulation algorithms in conjunction 
with a straightforward discretization scheme. 

Due to the promising results obtained by f inite=element algorithms in other 
fields of continuum mechanics and from initially promising applications to 
acoustics, finite elements were the medium chosen for this study. 

The great majority of development of the finite=element method has taken 
place within the context of structural mechanics. Concerning the relative 
advantage of various types of elements, it is informative to note the 
comments of a specialist in this field: Zienkiewicz (Ref. 18,^ page 104) who 
states that, "The question may well be asked as to whether any economic or 
other advantage is gained by increasing the complexity of an element. The 
answer, here is not any easy one, although it can be stated as a general rule 
that as the order of an element increases so the total number of unknowns 
in a problem can be reduced for a given accuracy of representation. 

Economic advantage requires, however, a reduction of total computations 
and data preparation effort and this does not follow automatically for a 
reduced "number of total variables as, though equation solving times may be 
reduced, the time required for element formulation increases. In general, 
the optimum element may have to be determined from case to case. 11 

In the light of these comments and for the sake of simplicity, a linear 
rectangular element from Zienkiewicz's Serendipity family (Ref. 18, page 107) 
was chosen as the basis for the analysis. Figure 2, contains a description 
of the nodal numbering system for the element, and shows orientation of 
local ( »2 , |)coordinates and global ( r ; Z) coordinates. (Since we assumed 
a harmonic angular ( 0 ) dependence by equation (16), only a two-dimensional 
discretization is required. ) 




A uniform discretization mesh for a duct of constant cross-sectional area 
is shown in Figure 3(a). The matrix map generated for this discretization 
by a nodal numbering system such as that shown in the figure may be 
developed as follows. Within each of the linear elements, interactions 
among nodal variables will only occur between those variables at the corners 
of the rectangle. Thus in element (2, 4), four sets of linear equations will 
link the four sets of nodal variables located at nodes 14, 15, 24, and 25. 
Consider the equation set representing variables at node 25. On assembly 
of the global matrix, this set will contain contributions from four adjacent 
elements namely (2, 4), (2, 5), (3, 4), (3, 5). That is, in the global matrix 
the variables at node 25 will be linked explicitly with variables at nodes 14, 
15, 16, 24, 25, 26, 34, 35, 36. All other coefficients in the submatrix 
representing the equation set for variables at node 25 will be zero. For the 
case of " £ " parameters per node, the matrix map generated by linear 
elements and the rectangular discretization numbered as shown in Figure 
3(a), is given in Figure 4. 

The principal generalized characteristics evident from the map shown in 
Figure 4 may be summarized as follows: 

• Global matrix consists of three diagonal bands. 

• Width of each band is "3 £ ". 

• Global matrix is block tridiagonal. 

• Major blocks are themselves block tridiagonal. 

• Order of major blocks is " nf " where 11 n " is equal to the 
number of nodes per row in the serial numbering direction. 

• Order of minor blocks is " £ ". 

Although the characteristics of a matrix conforming to a generalized map 

such as that described above are extremely convenient for minimizing 

storage requirements, and computational effort during assembly and 
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solution, consider whether much practical importance arises from this fact. 

The uniform rectangular discretization scheme from which this matrix was 
generated is hardly applicable to realistic problems. As shown in Figures 
3(b) and 3(c), however, the identical matrix map is also generated by a 
nonuniform rectangular discretization, as well as by arbitrary two-dimensional 
mappings of the discretization grid'T provided no fold-over occurs. 

Solution of equations (17) through (20) may be greatly facilitated by limiting 
the analysis to use of this discretization scheme. The specific map for this 
scheme may thus be incorporated implicitly into the computer program which 
enables us to disregard the general but somewhat ponderous sparse matrix 
pointer technology. All computation and storage is thus performed implicitly 
only with nonzero values. In addition, this discretization scheme yields a 
block tridiagonal matrix of block tridiagonal matrices. Special numerical 
techniques, which are discussed later (section 2.8), exist for solution of 
matrix equations involving tridiagonal matrices. 
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Element Derivation 


2. 3 

Consider the region of space enclosed by the element shown in Figure 2. 
Suppose that within this region the acoustic variables u.,1 as well as 
the mean fluid parameters T.tt.V.e of equations (17) through (20) 
are defined by relations of the following type 

1 U, 

£ 

where W, # U a f Uj J are the values of " U 11 at the four node points at the 
numbered corners of the element. 


= (*)!«} 


( 21 ) 


Define a set of local coordinates ( H ) within the region such that the coordi- 
nate of the nodes ( fji / ^ £ ) are given by (-1, -1), (+1, -1), ( + 1, +1), (-1, +1) 
for i = 1 through 4 respectively. The N.'s of equation (21) for the linear 
Serendipity element are then defined by 


^ ( , + ? ?‘) I' + ffr) 


( 22 ) 


For the undeformed element shown in Figure 2, the transformation from local 
to global coordinates is simply 

? , u^nl 

b (23) 

1 = 

CL 

where ( T c t Z c ) are the global coordinates of the center of the element. 
However, as illustrated in Figure 3, we require to map this simple shape to 
conform to complicated boundaries so that it then appears as shown in 
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Figure 5. Suppose then, that instead of equation (23) we have the transfor- 
mation ^ . 

ro 


r 
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(N, , N. , Wj , 



(K. W,, 



= Mir} 

= 


(24) 



Note that in comparing these relationships with those of equation (21) two 
conditions are satisfied: 


• The same points (viz. the node points) define the geometry 
and the finite element analysis points. 

• The shape functions (viz. the N^'s) defining the global 
coordinates and the variables within the element are the same. 


Thus, the mapped elements are isoparametric. 


In order to transform equations (17) through (20) into the local coordinate 

d d 

system, we need to derive expressions for , and - - 

dr dz 

system 


in the local 


3 Nj _ 3 Ni 3 r + 3 Nc dz 

3 1 ^ 9r 3 j 9z 3^ 


l 


3 A h — 3 A/i 3r t -f 3 A/t 3 z 

3 | d r 3 ^ 3z3| 







where 



is the Jacobian Matrix. 


(25) 


We may thus define the matrix 



(26) 


Suppose equations (21), (24) and (26) are substituted into equations (17) 
through (20). Using the Galerkin version of the method of weighted 
residuals, we wish to integrate over the element »' weighting each equation 
in turn by each of the N.'s to obtain the following four matrix equations in 
the nodal variables. It is convenient in integrating these equations to do so 
in the local, rather than in the global, coordinate system because this 
greatly simplifies the process due to the convenient limits of integration. 
Effectively, integration is then over the simple shape of Figure 2, rather 
than the complicated, distorted shape of Figure 5. In order to perform this 
operation, it is necessary to determine an expression for an element of area 
(dA) in the local system. It is shown in Appendix 1, that "dA" is given by 




d A = | J| d n ^ J 


(27) 


where ) X| is the determinant of the Jacobian matrix. Equations (17) 
through (20) thus become: 
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2.4 


Algebraic Integration for Ducts of Constant Cross-Sectional Area 


If we assume that we are dealing only with ducts of constant cross-section 
and that 

• the mean flow is incompressible, 

• there is no mean flow in the radial direction, 

• the mean axial flow has no gradient in the axial direction. 


Equations (17) through (20) how become 

= 0 


a«. / df> _ 


ku + v t N — — + 

dz B d2 

kv t M !*: + -L-it = o 
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From equation (23), we have that 
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whence it follows that 
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(36) 
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(39) 


( 40 ) 



Supposing a linear dependence for M within each element, 

M = M c ■+ M'(r- r c ) 

= M c •* M' ij 

We may integrate over the element using the Galerkin procedure as in 
section 2. 3 (note that now IJI , the determinant of the Jacobian, is unity) 
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where (i = 1 to 4) yielding a total of sixteen equations. Writing U.jVjW' and jy 
in terms of the nodal coordinates as in equation (21) 


u.= (.N) 5.U.} 

v = (v) {/} 
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equations (41) through (44) may readily be integrated. 


Before performing these integrations, however, let us consider how 
conditions at the element boundaries may affect the integration. If 
equations (32) through (35) contained second-order terms it would be 
necessary to integrate these by parts to reduce them to first order since 
second-order terms require continuity of first-order terms over a 
boundary (Ref. 18, page 42). However, using a linear element with a 
linear system of equations has advantages in this respect since we only 
need to satisfy the requirement of continuity of the variables themselves 
across boundaries. Also, along the outer boundary, it is sufficient to 
constrain the appropriate variables to conform to the boundary conditions 
without the necessity of evaluating the specific boundary integral. 


Equations (41) through (44) become 
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m - ipi ■ m ■ in - I>] are defined in appendix 2 . 
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2. f> Numerical Integration for Arbitrary Axisymmetric Ducts 

Numerical integration of equations (28) through (31) by a standard inte- 
gration formula is a trivial procedure in concept. However, because of the 
necessity to integrate a large number of elements in this manner, consider 
the equations more carefully and optimize the process. 


Since none of the variables is an explicit function of " Q ", the integration 
with respect to " 6 " results simply in a multiplicative constant of 2 
which may be ignored. 


It may be seen that many coefficients of the matrix equations are multiplied 


by sc 

:alar quantities. 

These are ( 

defined below: 

R = 

W lr } 

8 

= (*) £J ? 

c * 

OW {V} 

P 

= M {« } 

B = 


F 

= 

Sc = 

Miui 

a 

= («^ z ) {U } 

X - 

Nil} 

L 

- (^r)(c) 

K = 

M £< i 

* 

r*>~\ 

WJ 

It 

M = 

W J 




(50) 


where the quantities jf) through G? are functions of ^ and . 

Substituting in equations (28) to (31) we get 
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The Gaussian quadrature procedure that is used employs unevenly spaced 

abscissa points ( The corresponding weights to be applied at these 
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Tables of abscissae and weights for different values of n are given in 

G 

and are reproduced in Appendix 3. 

The logic diagram for evaluation of equations (51) through (54) using 

t' 

equation (55) is shown in Figure 6. It may be seen that it is possible to 
retain several potentially time-consuming steps outside the innermost loops 
of the quadrature process. 

2. 6 Global Matrix Assembly, Packing Technique, and Insertion of 
Boundary Conditions 

Assembly of a global matrix equation, that is, a matrix equation represent- 
ing the entire system, from a set of finite element matrix equations is basic 
procedure in the finite element technique, e. g. , reference 18, page 13, 
and is not described in detail. From Figure 3(a), it may be seen, however, 
that there is a one-to-one correspondence between the global serial number- 
ing system of nodes and the element numbering system together with its 
local node numbering as shown in Figure 2. Appropriate shifting of rows 
and columns are all that is required to add the local element matrix directly 
into the global matrix. 

This procedure is somewhat complicated by the fact that when a global 
matrix possesses.® map^as that shown'in Figuf^ejAyeitlis undesirable' 
to store anything but the three diagonal bands of nonzero coefficients. The 
most convenient way to "pack" a matrix of this form is shown to the right 
of the full or expanded matrix in Figure 4, where the diagonal bands are 
arranged in vertical columns. The algorithm to yield a column number in 
the packed matrix from a row and column number in the expanded matrix is 
listed in Appendix 4. 

Along the outer boundaries of the discretised region, two kinds of boundary 
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conditions may exist, noise source or admittance. The boundary condition 
at the radiation plane may be forced to conform to an admittance boundary 
condition as shown in Appendix 5. 

Provision is made for insertion of the two kinds of boundary condition into 
the assembled global matrix equation by imposing relations on the nodal 
variables along each boundary. 

The noise source boundary condition consists simply of setting all nodal 
pressure coordinates equal to a specified constant value (implying by equation 
(16) a specified harmonic amplitude). For example, see Appendix 6. 

The admittance boundary condition involves specifying a linear relation- 
ship between nodal pressure coordinates along the boundary and the 
normal component of acoustic velocity. For example, see Appendix 7. 

The examples in appendices 6 and 7 assume that, providing the global 
matrix is stored in packed form, it may be held entirely within direct access 
memory. This is not the case however, because by definition (see section 2. 2) 
we will be dealing with extremely large problems employing over one hundred 
thousand degrees of freedom. 

During assembly of the global matrix, however, it is unnecessary to have 
the entire matrix in direct access memory. Figures 3a and 4 show that contri- 
butions to the equations for variables at global node number 25 (for example) 
may only come from elements (2, 4), (2, 5), (3, 4), and (3, 5). Thus, it is 
only necessary to hold variables corresponding to two element rows (or three 
nodal rows) in direct access storage at any particular time. This is also the 
case for application of boundary conditions, since in Appendix 7, equation (1), 



for example, > x t , and X m are all variables attributable to a single 
node and thus all lie within the diagonal submatrix. . Suppose that their equiva- 
lents are specifically attributable to node 25 in Figure 4. The only rows and 
columns affected by application of the boundary conditions thus lie within the 
bounds of those attributable to node 25. Thus, if matrix block rows from 
number 1 1 to number 30 are in direct access memory, the boundary condition 
may be applied to these only with the result being identical to an application 
to the entire matrix. 

2. 7 Solution of Matrix Equation 

The matrices generated by the Galerkin formulation of sections 2.4 and 2. 5 
are in general unsymmetric and nonpositive definite in addition to being 
complex. Due to the discretization scheme described in section 2. 2, they are 
block tridiagonal. 

The JLU -decomposition of a block tridiagonal matrix is readily obtainable 
(Ref. 19, pages 166 and 173). The L,U -decomposition takes the form 
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where m is the identity matrix and where the two matrices on the right- 
hand side of the equation are the lower [ijand upper [ul triangular factors, 
respectively. The component matrices of ft I and [u] are given by direct 
substitution as 
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which may be solved in turn in the usual manner by forward and backward 
substitution using the relations 
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The process described in equations (56) through (67) may be used to solve 
the global matrix equation whose derivation and assembly is described in 
this report. We described in section 2. 6 however, that assembly of the 
global matrix equation may be carried out efficiently by keeping only three 
nodal rows of coefficients in direct access memory at any one time. 


In terms of the required size of direct access memory, three nodal rows of 
coefficients reduces to nine packed blocks and three block vectors (e. g. , see 
Figure 4). The respective orders of these matrix elements are: 
packed blocks fy\ H % 3^ 

block vectors YY\ i X 1 

where m = number nodes in serial numbering direction 
__ £ - number of parameters at a node'.'; 

• ,Thus, the total storage requirement for assembly i s ^ 2. 7 W -c + 3/7? C J , 


For solution of the matrix equation, however, an efficient minimum require- 
ment is three packed blocks, three block vectors, and one expanded block, 
giving a total storage of ( ^ W? 0 ~b 3m i •+ The necessity for the 

storage-consuming expanded block is clear from equation (58) where it may 
be seen that [Pkl is a full matrix for all ( k - 1, PI ) • 



Careful inspection of equations (57) through (67) shows that a number of 
specialized matrix handling algorithms need to be used in order to limit the 
requirement on the number of expanded blocks to one. These include the 
following: 

• Transposition of either a packed or expanded matrix, and 
superimposing the transposed matrix upon the original storage 
area. 

• Multiplication of an expanded matrix by a packed matrix, and 
superimposing the result over the storage area of the original 
expanded matrix. 

• Addition and subtraction of expanded and packed matrices, and 
superimposing the result over the storage area of the original 
expanded matrix. 

• Multiplication of a packed matrix by a vector without expanding 
the packed matrix. 

Once the necessity for these algorithms is recognized, however, their 
preparation is trivial. Their use in the implementation of equations (57) 
through (67) together with a map of direct access storage requirements is 

i ~ > 

shown in Figure > 7, 

A 

<** / 

2. 8 Acoustic Attenuation 

The principal final result of an analysis of sound propagation in an engine 
duct is a value for its acoustic attenuation. 

Clearly from equation (16), any evaluation of the total attenuation for sources 
emitting sound over a broad frequency band, and with an arbitrary spatial 
distribution of amplitude and phase over the source plane will need to be 
carried out over the entire double Fourier transform implicit in this relation- 
ship. This implies practically that a value of cO has to be specified for 
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Figure 7. Description of Block Tridiagonal Matrix Solution 
with Memory Map. 



































each time harmlonic", “‘and a value' of W for each-angular harmonic. 

Since the [idjM) harmonics are orthogonal, acoustic intensities computed 
from them independently have simply to be summed to produce total intensity. 


The choice of an expression for acoustic energy is not immediately clear. 
For irrotational, uniform -entropy flow (in the Z -direction) , Cantrell and 
Hart^^ , and Morfey^^ have shown that 



( 68 ) 


where I is the axial intensity, and * denotes complex conjugate. The above 
z 

expressions are accurate to second order in fluctuating quantities. However, 

both Cantrell and Hart, and Morfey obtained their expressions from the time- 

( 22 ) 

averaged energy equation. Eversman used, as his starting point, the 
exact energy equation for a nonviscous, nonheat-conducting flow and obtained 
for the axial intensity: 


<I> 



where V - radial component of acoustic velocity 

W - tangential component of acoustic velocity 
Both results were programmed, with somewhat more consistent answers 
being obtained from the second expression. (In some instances, use of either 
expression results in a net increase in flux at the radiation plane over the 
flux at the source plane.) 
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In order to obtain axial acoustic flux over source F and radiation F 

s r 

planes, a simple linear integration scheme described in Appendix 8 was 
used. Duct attenuation is thus given by 
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3.0 DESCRIPTION OF RESULTS 


3. 1 Comparison of Model with Previous Analyses 

Comparison of the model with previous work is complicated by the fact that 
these comparisons can only be achieved using a small proportion of the 
model's total capability. These comparisons do have the advantage, however, 
of providing some confidence in the convergence characteristics and overall 
accuracy of the model. 

A uniform cylindrical duct of unit radial and axial dimensions with an outer 

wall impedance of pc is shown in Figure 8. Duct attenuation was calculated 

at a wavenumber of unity in the absence of mean flow at a circumferential 

harmonic number of zero for a series of element assemblies. The solutions 

obtained from the model may be seen to converge on a suliton for the same 

( 15 ) 

problem computed from the work of Zorumski 

In order to verify the model's accuracy over a wide range of parameters, 
and to evaluate its feasibility for optimization studies, two optimization 
schemes were programmed. The first was a two-variable optimization 
to determine real and imaginary parts of a uniform liner admittance yield- 
ing maximum acoustic attenuation. This optimization was performed at 
two frequencies and the results are compared in Table I with those derived 
by Lester and Posey and Quinn . The results from the three 
analyses are substantially the same. A six- variable optimization for a 
three-section liner was then carried out with no significantly different results 
from those of Quinn. These results, together with those of Quinn, are also 
given in Table I. 

The optimization scheme used in the two determinations above was the 
Davidon Fletcher Powell technique adapted to numerical gradient compu- 
tations. The fact that reasonable agreement with previous analyses was 
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Figure 8. Comparison of Acoustic Attenuation vs Mesh Density with 
Exact Solution Computed from Work of Zorumski^, 






TABLE I 

Comparison of Optimum Acoustic Liner Determinations 


Duct radius = lm Duct length = 2m c = 344m/ sec 

Circumferential wavenumber = 1, no flow 



548. 2 
344.4 

548. 2 


Liner Admittance ( x 1/pc ) 




Lester 


Quinn 


Duct -Attenuation ,(dB) 


t _ „4.__(23) 



Le ster 


Quinn 


. 21+ . 45i 

. 25+ . 38i 

. 27+ . 55i 

. 29+ . 56i 


0+ . 40i 
. 16+ . 57i 
. 60+ . 59i 


UNIFORM LINERj 

4. 7 
7. 3 


THREE -SECTION LINER 


. 003+ .49i 
. 31+ . 46i 
. 71+ . 67i 



TABLE II 

Comparison of Duct Attenuations in Presence of Plug Flow 


Duct radius = lm 


Duct length = 2m 


Circumferential Liner 
Hz Mach No. Wavenumber Admittance 

( x 1/pc ) 


c = 344. 4m/ sec 


Duct 

Attenuation (dB) 


Current 

Analysis 




Zorumski 
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obtained is encouraging because evidently a sufficiently low "function 
noise" was present in spite of a moderately coarse analysis of 12 x 30 
elements. 

A conical horn was selected to verify the program for nonuniform geometry. 

A solution based on a small angle approximation at low frequencies is given 
(25) 

in Morse . This solution assumes a spherical wave emanating from an 
equivalent point source at the produced apex of the cone. The amplitude 
of the pressure is thus inversely proportional to the distance from this 
hypothetical source. Figure/9 . 'shows the pressure variation with axial 
distance along a radius from the hypothetical source selected arbitrarily 
at 12. 9° to the axis. Results from the current model (shown by a solid 
line) show little deviation from results produced by Morse's analysis 
(open circles). 

A. limited amount of published data is available for ducts containing flow 

that is suitable for comparison with the current model. In the absence 

(15) 

of an alternative comparison, Zorumski's analysis was used. This 
comparison is less than ideal in two ways. The current analysis is a 
shear-flow model, that is, the mean flow is required to go to zero at 
duct walls whereas Zorumski's is a plug-flow model. Further, the current 
model has a p‘c termination impedance, whereas the version of Zorumski's 
model used, had a zero-reflection termination boundary condition. 

Reasonable comparison is demonstrated from the results presented in 
Table II. For these analyses, plug flow was simulated as closely as 
possible in the current model by a very small shear layer. Several 
combinations of frequency, Mach number, circumferential harmonic 
wavenumber, and liner admittance are shown in Table II. No significant 
deviations were observed in computed duct attenuation. 
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3. 2 Illustjra tions of Potential Current Model in Duct Optimization 
The current model presents the opportunity to explore novel concepts in 
aero-engine duct design for optimum acoustic characteristics. This 
section contains illustrations of three such concepts and describes their 
relative success. 

3. 2. 1 Quadratic Liner Variation 

As an alternative to a sectioned liner, a uniformly varying liner 4s 
suggested. To compare directly with the three- section liner whose 
results are given in Table I, a quadratic liner variation was programmed. 
As a viable alternative in optimization studies to a three- section liner, a 
quadratic liner possesses the same number of variables. Its admittance 
3 is given by 


B(z) = a + a_ z + a_z 
o 1 2 


where the a's are complex. 


For the duct similar to that described in Table I, optimization coefficients 

) ■ t-Q 

yielding an admittance function given in ^Figure 10 was derived. Duct 

' \ 'Ci C 

attenuation was 5.9 dB which compares closely with^he value of '6. 2 dB 
obtained for the three-section liner described in Table I. 


3. 2. 2 Optimization of a Convergent-Divergent Duct 

In addition to acoustic-liner parameters, the physical shape of the duct 
itself may be optimized for maximum acoustic attenuation. Suppose, for 
example, a constrictipn defined by a cosine function is imposed on the outer 
^ wall. The duct radius is :given;by 


R = j^l - q(l - cos 2TTz)j- R q 
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where q is the constriction oarameter and R is the duct radius at z = 0. 

o 

Suppose further that we decide to place a lower limit on the constriction 
size such that cross-sectional area at the throat should be no less than 
one-half the cross-sectional area at the source plane. Thus q should be 
in the range 

0 < q < * 15 

This may be achieved by the transformation 

q = . 075(1 + cos Q) 

A three- variable optimization was carried out with the real and imaginary 
parts of a uniform admittance as variables in addition to the duct constric- 
tion parameter Q defined above. 

The admittance was found to converge rapidly to the same optimum value 
which was previously determined for a similar duct with no constriction 
at the same frequency, and circumferential harmonic number. The 
constriction parameter, as might be expected, converged to the limiting 
value of 


q = . 15 

While this exercise may seem somewhat trivial, it serves to demonstrate 
the possibilities of duct optimization using the model. 

3, 2. 3 Optimization with a Center Body 

The shear layer produced, as the mean flow veloCity. r goes to ze ro ^ori 3 * 5 - the-" 

walls of an aero-engine inlet duct, causes-acous'.tic eiTergy t-o be-refract^d 
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towards the duct axis is well known. This effect was verified using the 
model and is depicted graphically in Figure ll'T 

An unfortunate consequence of this effect is the decreased effectiveness 
of acoustic liners placed on the outer duct walls. To use an acoustically 

lined centerbody, jsuch as that. shown in Figure 12, shaped so as' to 

\ ^ 

minimize its o.wn refraction characteristics, ,.w.ould:,app:ea.r profitable. 

A three- variable optimization was carried out on a duct of this type using 
a cosine centerbody function similar to the function described in the pre- 
vious section for a constriction. The outer duct wall was uniform 
(radius = 1 m), and the maximum radius of the centerbody was allowed 
to vary between zero and . 6 of the outer wall radius. Duct length was 
2 m. The other two variables in the optimization were the real and 
imaginary parts of the liner admittance on the centerbody and the outer 
wall. 

A mean flow model based on streamlines, calculated on the assumption 
of incompressible flow, was used. Free stream velocity at the source 
plane was Mach . 2, with a one-half sine rolloff profile. In the absence 
of better information, this profile was assumed to have a constant normal- 
ized shape which was independent of axial station. 

At a frequency of 548. 2 Hz, the centerbody radius increased to an optimum 
at its maximum constrained value of . 6 m. The corresponding liner ad- 
mittance, yielding an optimum duct attenuation of 25. 8 dB, was 
(. 62 + . 38i). 

This attenuation appears remarkably large when compared with the corres- 
ponding values in Table I. A partial explanation, however, is undoubtably 
that the dimension between the centerbody and outer duct wall is of the 
same order as the wavelength of the first mode. 
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Figure 12. Constant Cross-Sectional Area 
Inlet Duct with Center Body. 
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3. 3 Computing Time 

The finite element model described in this report was specifically designed 
to operate making most efficient use of computational resources. Thus, 
perhaps the most important result of the study is the degree to which 
this objective is fulfilled. 

Figure 13 shows typical central processor times that were achieved 
from the different variations of the model. We may draw the following 
conclusions from the results: 

• For a particular matrix bandwidth, central processor time 
is roughly proportional to the matrix order. 

• Numerical integration increases total computing time by roughly 
a factor of 1. 3 over the algebraic form. 

In comparison with NASTRAN (NASA Structural Analysis), the general 
purpose structural finite element program, the current model shows a 
speed improvement for similar matrices of approximately a factor of 12. 
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CONCLUSIONS 


4 . 0 

The objective of this work was to develop a finite element scheme for 
the acoustics of aero-engine ducts. The principal difficulty of a direct 
numerical solution is containing the problem within existing computers 
and attaining a solution in a reasonable amount of computing time. The 
mathematical formulation, logical flow sequences, matrix manipulation 
techniques and a solution scheme to attain the stated objective are described 
in this report. 

Testing and verifying the model is complicated by the fact that available 
alternative analyses possess only a small portion of its total capability. 

With this restriction in mind, testing of the model was carried out in such 
a manner in order to compare results with several alternative analyses. 
Results of the testing program are satisfactory with the model showing good 
convergence characteristics over a wide range of conditions. 

Specific parameters such as spatial acoustic pressure variations and duct 
attenuation were compared with alternative analyses. General trends such 
as refraction of acoustic energy by a mean flow shear layer were verified. 

An aspect of the model likely to yield interesting information on acoustic 
design of aero-engine ducts is its ability to be used in conjuction with 
optimization studies where parameters such as duct shape, size and 
shape of centerbodies, and multisection or continuous acoustic liners may 
be varied. Illustrations of the use of the model in this context were presented 
together with results showing expected and unexpected trends. 
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Testing of the model in the presence of compressible mean flow has not 
yet been accomplished. This capability is present in the model but requires 
linking with a separate analysis to provide mean flow data. 
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Appendix 1 


Area Element Coordinate Transformation 


Suppose a vector R is given in the global ^Z t r) system as 
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Appendix~2 
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. 2- 1 

Definition of terms used -in equations (46) through (49) ^ 
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wlrere for example, is the 'th component'of the 4 x 4 matrix [A] . 





Abscissae and Weight Coefficients of the Gaussian Quadrature Formula 


Jw dx 

-I 

Hj f(*j) 






j*/ 

i " 










H 



n 

= 

2 



0- 57735 

02691 




1-00000 

00000 



n 

= 

3 



0* 77459 

66692 




0- 55555 

55555 

0- 00000 

00000 




0-88888 

88888 



n 

= 

4 



0-86113 

63115 




0-34785 

48451 

0-33998 

10435 




0- 65214 

51548 



n 

= 

5 



0-90617 

98459 




0-23692 

68850 

0-53846 

93101 

‘ 



0*47862 

87604 

0-00000 

00000 




0- 56888 

88888 



n 

= 

6 



0-93246 

95142 




0- 17132 

44923 

0-66120 

93864 




0-36076 

15730 

0-23861 

91860 




0-46791 

39345 



n 

= 

7 



0-94910 

79123 




0- 12948 

49661 

0-74153 

11855 




0-27970 

53914 

0-40584 

51513 




0-38183 

00505 

0-00000 

00000 




0-41795 

91836 



n 

= 

8 



0-96028 

98564 




0- 10122 

85362 

0-79666 

64774 




0-22238 

10344 

0-52553 

24099 




0-31370 

66458 

0-18343 

46424 




0-36268 

37833 



n 

= 

9 



0-96816 

02395 




0-08127 

43883 

0-83603 

11073 




0-18064 

81606 

0-61337 

14327 




0-26061 

06964 

0-32425 

34234 




0-31234 

70770 

0-00000 

00000 




0-33023 

93550 



n 

= 

10 



0-97390 

65285 




0-06667 

13443 

0-86506 

33666 




0-14945 

13491 

0-67940 

95682 




0-21908 

63625 

0-43339 53941 




0-26926 

67193 

0-14887 

43389 




0-29552 

42247 


-59 



Appendix 4 

Listing of Subroutine PINDS 
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Appendix 5 

Termination Boundary Condition 

Assume that in the radiation plane U varies slowly over an element and 
that a plane wave solution holds, i. e. , 

-ikz 

U = UL 0 £ 

-ikz 

W- « 

where k = M. - 

a c + a 

The axial momentum equation (equation (6) was 

Jc<Ju + (fr+f7)iU + 77i£ 

dr dr 

+(fu + ft I)iu +TR2U + 2± 

9z. 

assuming that 9U , 1/ , 5^ = 

5z. 


this reduces to 


SilPu. + J u 2 jl + iIl = ° 

?Z dz 

Substituting (1) and (2) into (4) gives 
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Equation (5) is an admittance condition, and is the approximate boundary 
condition applied at the radiation plane. 
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Appendix 6 

Insertion of a Noise Source Boundary Condition 


Consider the following system of equations!. 

'm N = { 8 } 
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/ Appendix 7 


- 


Insertion of an Admittance . Boundary Condition 


7-1 


Consider the- following system of equations: - r -v 

WW= {8} 


•t ^ 

i; e.- 
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The boundary condition P-f- , where p is the wall admittance and Vq 
is the normal component of the acoustic velocity, may be rewritten as 


V = C, + Cj. K 

, = t 

where ' cos 

Cj. = t°.n{e) 

and 0 is the angle between the wall and the "Z-axis (see Figure 7 "q). 


Suppose 



then 




C/ x m + C 4 X* 



( 4 ) 
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7-3 


5 • A- - f/ 

Substituting (4) in (1) gives 

^//X, v- — + {CL lk +c x a.,i))(k o.*t -„ + (a lni +c,4, e )x^„+a ln x n = l, 
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( 


Clearly, we may have chosen to eliminate either (A or p instead of V 
in the process above. The choice of which parameter to eliminate depends on 
the numerical values of C, and . If in equations (3), for example, 9 

tended towards jT , unreliable results will be obtained from equations (5). 
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Appendix 8 

Linear Numerical Flux Integration 


In order to complete acoustic flux it is necessary to evaluate the intensity 
distribution over a plane. This may be done numerically in the following 
manner (assuming ( -r r, ) is smallf : 
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